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ABSTRACT 

Starting data are provided for internal characteristic programs 
pertaining to axially symmetric rocket nozzles. Several approximate 
methods are examined before a complete analysis is attempted. 

Using perfect frozen gas considerations, the potential function in 
the region of interest is approximated by a double power series in the 
space variables. The potential equation of motion, with the inviscid 
boundary condition, produces non-linear simultaneous equations. 

The non-linear equations are handled uniquely, and the results 
are utilized to describe the flow field. Different methods of checking 
the validity of the results are applied, and comparisons are made 
with other analyses. 

The remaining difficulties in the development of a production 
program that can handle arbitrary minimum section geometry and 
slightly varying mass flows are discussed. 
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SYMBOL 

a 


a* 


a 


i, 


n 


A ! 


2j 


A(i,2j) 

A i 

Bi 


b i 
Di 
f 2i 

g(y) 


1. J 


G(m, p) 


G(i, p) 


DEFINITION OF SYMBOLS 

DEFINITION 
local speed of sound 

local speed of sound where the velocity is sonic 

coefficients in a series expansion of g(y), or 
constant in Method I approximate analysis 

coefficients in the polynomial representation of 
G i.P 

coefficients of the double power series expansion of 
the potential function, 

same as A ^ 2j 

coefficient used in approximate analysis II 

coefficients of the series expression for the bound- 
ary curve, or coefficient used in approximate 
analysis - Method II 

coefficient used in approximate analyses in Method I 
for f 2 i while in method II for a^ 

the "unknown" of the i tb boundary condition 
equation, i.e., A i+ i , 0 

coefficients of the $ series in the approximate analyses 
x = g(y) defines the surface of the sonic line inMethod II 
indices 

denotes any of the boundary condition equations (for 
m <M) containing p+ 1 columns 

4.-L 

denotes the i L boundary equation, and contains 
p+ 1 columns 
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DEFINITION OF SYMBOLS (CONTINUED) 


SYMBOL DEFINITION 

G(i,oo) a function equal to zero, denotes the i ** 1 boundary 

equation and contains all of the columns 

m mass flow (dimensional) 

mass flow (non-dimensionalized) 

denotes the number of D. unknowns considered; 
it is one less than the number of rows considered 


& 

M 


M* 


Mach number referred to the speed of sound a* 


m 


index of the boundary condition equations 


n 


index 


P 

s 

u, V 
U, V 


X, y 

a 

e 


P 

e 

y 


integer, one less than the total number o£ columns 
considered 

subscript denoting boundary 

perturbation velocities in the x and y directions 
respectively, non-dimensionalized on a* 

total velocities in x and y directions respectively, 
dimensional 

minimum throat radius 

axial and radial directions in circular cylindrical 
coordinates 

constant, equal to axial velocity gradient 

distance from minimum section to the sonic line 
along the x axis 

throat radius of curvature at the minimum section 
longitudinal coordinate, £= x-g(y) 
ratio of specific heats 
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DEFINITION OF SYMBOLS (Concluded) 



i ► y ~ l 

equal to 

perturbation velocity potential 

partial derivatives of <J? , equal to u and v, respectively 
flow angle 

zero or one function 
isentropic stagnation density 
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SUMMARY 

1 1 

This study will ultimately produce a working computer program 
to provide the initial data for the method of characteristics. The 
program must be more exact than those analytical methods now 
available. 



Since the flow becomes supersonic in the vicinity of the throat of 
the nozzle, analysis of the transonic flow in the throat region is of 
interest. The physical problem was approximated by a steady-state 
axisymmetric system (co-ordinates x and y) with zero radial velocity 
along the axis. In this region, the fluid is considered to be a homogen- 
eous, perfect gas in frozen equilibrium with no viscous forces. The 
system is also considered irrotational and isentropic with negligible 
energy losses through the nozzle surface. In the region of the throat, 
the nozzle surface is expressed as a power series in terms of the 
axial coordinate and is not limited to boundary surfaces of constant 
radius of curvature for the minimum section. 


This method attempts to solve the full, non-linear potential 
equation of motion. A double power series in both x and y is assumed 
for the potential function. The original equation is satisfied by re- 
currence relationships for the general coefficients that are evaluated 
in terms of the coefficients of the velocity distribution along the axis. 
The remaining coefficients are determined by satisfying the inviscid 
boundary condition along the nozzle contour, where the contour i 
described by a power series in x. 



The problem is that the analysis produces an infinite set of highly 
non-linear simultaneous boundary condition equations. The number of 
equations used was determined by the size of the region of convergence 
of the flow field. The minimum number was solved by subsequent 
iteration. The iterative method deliberately does not satisfy each 
equation. The error of the equation is used as a bias so that the effect 
of the iteration of the rest of the system will channel the forced error 
toward a correct solution. This method should be usable in other 
areas where non-linear simultaneous equations occur. 

A library of solutions is needed to pick initial input values close 
to the actual solutions. The results are presented for constant radius 
of curvatures of 5 to 1. 25 times the minimum section radius. 

This analysis may be the first to allow small variation in mass 
flow in a given nozzle. The variation appears as a shift in the loca- 
tion of the point where theMachone line intersects the axis. 

The solution is used to compute input data for characteristics 
and to check the accuracy of the program. Although this method is 
a frozen flow analysis, it is used in conjunction with equilibrium 
programs by stipulating the value of the ratio of the specific heat from 
an equilibrium combustion and expansion program. The method of 
characteristics was used to confirm the validity of this method. 

The obstacle to general use of this method is the amount of com- 
puter running time, but the program should prove to be valuable in 
nozzle design. 
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INTRODUCTION 


The method of characteristics, as applied to gas dynamic problems, 
is restricted to certain regions in the flow field. These regions are 
those in which the flow is not only supersonic but also purely hyper- 
bolic in the mathematical sense. Practically, it is restricted to 
those regions in which the local Mach angle is not large enough that 
the resulting characteristic mesh calls for excessive computer time. 
Therefore, it is necessary to generate, by some other analysis, 
certain calculated flow properties over a suitable surface from which 
the method of characteristics can be started. 

Ultimately this study will produce a working computer program 
that will provide the starting data for the method of characteristics as 
applied to the flow field internal to a rocket exhaust nozzle. Since the 
flow becomes supersonic in the vicinity of the throat, the analysis of 
the transonic flow in the throat region is of interest. There are 
several approximate methods for treating this particular flow; however, 
none of these have an accuracy obtainable by characteristics (1)*. 

Some of these methods are discussed to show their salient features, 
regions of applicability, and limitations. 

Also discussed in this report are the method developed to calcu- 
late the flow field to within an arbitrary degree of accuracy, the 
computer program utilizing that method, and a comparison of various 
methods. 

This study was performed by R. S. Mendelson of the Structures 
and Mechanics Department of Chrysler Corporation Space Division, 
Huntsville, Alabama Operations under contract NAS 8-4016, Task 
Order M-P& VE-PA-M-6, Task Assignment M-P& VE-PA-6-63. 


^Bracketed numbers that appear as supperscripts refer to the 
Bibliography. 
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ANALYSIS 


1 . GENERAL 

There are several methods, other than the two presented in this 
paper, which give approximate solutions for an internal transonic 
flow. Before discussing the new method, which is arbitrarily exact, 
some of these methods should be studied because they serve to check 
the exact method. 

In these analyses, the physical problem has been approximated by 
a steady state axisymmetric system. The fluid in the region investi- 
gated is assumed to be a homogeneous, inviscid, perfect gas of con- 
stant composition. Although there is energy lost to the nozzle wall, 
these losses are neglected, and the system is considered an irrota- 
tional, isentropic one. The boundary is required to be continuous 
and to have a continuous slope. The exact analysis also requires that 
the nozzle radius in the region of the throat be expressible as a power 
series in terms of the station variable x. 


Combining the equations of- continuity and momentum, and impos- 
ing the irrotationality condition results in the general equation of 
motion 


( a 2 - U 2 


3U 

a x 


- 2UV 


3U 
3 y 


+ (a 2 


- V 2 ) 


3 V 

a y 


+ a 


v=o< 2 > 

y 


(i.D 


where: a is the local speed of sound,. 

U and V are the total velocity components in the x and y 
directions, respectively. 

x and y are the axial and radial directions in circular, 
cylindrical coordinates. 

The speed of sound, where the Mach number is unity, is defined 
as a* and is only a function of stagnation conditions. The velocities 
U and V are nondimensionalized with respect to a* and written in the 
form of perturbation velocities, u and v, so 

-, = 1 + u (1. 2) 

a* 
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and 


— = v 
a* 


(1. 3) 


System dimensions, x and y, are nondimensionalized by the 
minimum nozzle radius y^ (FIG 1). 

- The equation of state and the isentropic relationship and equations 
1. 2 and 1. 3 are used to reduce equation 1. 1 to the form 


9 u 
d x 


( Y + 1) 


-r u (2 + u) - v 2 ] §f - 


ju(2+u) + r v 2 ] 

a T^7 + (yrr - r[.<z + »,♦.*]) = - 


4v 


(1 


where: 


7 is the ratio of specific heats, and 

y - 1 


r = 


7+1 


ay 

o (1.4) 
(1.5) 


The condition of irrotationality assures the existence of a velocity 
potential $ , so 


and 


$ = u 

x 


$ = v 

y 


(1. 6) 


where the subscript notation indicates partial differentiation with 
respect to x or y. 

By virtue of equation 1. 6, equation 1. 4 is reduced to the form 

J - (2+$ x )$ x - r ^yj $ xx + + j ) “ r $ x (2 + ®x) ~ ^y 2 | ®yy “ 

WTT) (1+ ®*> % ®*y + (tttT) - r K < 2+ ®*> + ®yl)4 I = 0 (1 ' 7) 

Equation 1. 7 is the basic equation from which all succeeding analyses 
originate. 

2. APPROXIMATE ANALYSES 

Deleting terms of the order $ x 2 and $y 2 as small compared to 
unity, equation 1. 7 can be simplified to 


(7+ U $X $XX + 2$y $xy - — J- (y$y) = 0 


( 2 . 1 ) 
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y and Y 


NOZZLE SURFACE 



FIGURE 1 LOCATION OF COORDINATE SYSTEMS 
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a. Sauer's Method (3) 


In Sauer's analysis, the origin of the coordinate system 
(FIG 1) is placed on the axis at a distance c (to be calculated later) 
downstream of the minimum section so that it coincides with the 
intersection of the axis and the critical (sonic) line. Hence, from 
equations 1. 2 and 1. 3 and the definition of equation 1. 6, it follows that 
$ x (0, 0) =$y (0, 0) = 0. Furthermore, it is argued that both $y and 
$ X y are small in the vicinity of the axis, y=0, and accordingly, the 
term 2$y <5 X y is deleted from equation 2. 1, resulting in 

( 7 + 1 ) -Syy- -y*Y = 0 ( 2 - 2 ) 

Sauer assumes the existence of a solution to equation 2. 2 in 
the form: 


$ 



2i 


(x) 



where the coefficients f 2i( x ) are functions of x only. 


(2.3) 


The derivatives of $ from equation 2. 3 are substituted into 
equation 2.2. The resulting recurrence relationship and the 
appr oximation 

f o ( x ) * ax 


where: a - constant, results in 

f o (x)= « if ] 


(2.4) 

(2.5) 


f2 (x) = 


(7+1) 

4 


or 2 x 


*4 (x) = 


(7+1) 2 a 3 
64 


with f 6 (x) = fg (x) = f n (x) = 0, n> 6 


The constant or is a first approximation to the longitudinal velocity 
gradient on the axis of symmetry. 


Substituting equation 2.5 into equation 2. 3 and differentiating gives 


and 


<5 x = u= a x + 


(7+1) 

4 


a 


$y = V= 


(7 + 




xy + 


(7+1) 2 a 3 y 3 
16 


( 2 . 6 ) 
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The two following boundary conditions are required to evaluate a 


and e : 


and 


at x = -£ and y = 1 
v = 0 


1 


3 v 


p (1 + u) 3 x 

where p is the radius of curvature at the minimum section. 


(2.7) 


Equation 2.7 insures that the streamline at the wall has the slope and 
curvature of the wall at the minimum section. Solution of equation 
2. 7 gives 


and 


y + 1 

e - — a 

8 


_1_ 

P 


(7+1) a z 


} ^ «*] 2 

Sauer assumes that in equation 2. 9 

Hr>’ 1 

and drops that term to get 


( 2 . 8 ) 


(2.9) 


a yjw TTI T 


( 2 . 10 ) 


But closer examination shows that terms of the order of a 2 have not 
been dropped any place else in the method. If the term is not dropped 
from equation 2.9 then equation 2. 10 becomes 


a = 


/ 2 - 

4 1 

/( 7+ 1) 

[4 p - lj 


( 2 . 11 ) 


Obviously, unless 4p>>l, the 7 * ■ ~ a 2 term cannot be legitimately 
dropped in equation 2.9. 
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The relationship 


(1+u) 2 + v 2 = l (2.12) 

defines the critical curve M# = 1, where M*=^(l+u) 2 + v 2 . 

Sauer used for a first approximation of the curve the condition 

u = 0. (2.13) 

Sauer's method was never intended to be an exact procedure for 
obtaining the solution of the transonic flow field. Consequently, 
Sauer's method should not be expected to produce results having 
accuracy comparable to the method of characteristics. In the dis- 
cussion of the other simplified analyses, various approximations that 
Sauer used are eliminated. 

b. Other Methods 


In addition to Sauer's approximation, there are other approximate 
analyses which give useful comparisons. One of the more notable 
analyses is that of Oswatitsch and Rothstein (2). It, along with that of 
Sauer, is discussed briefly by Shapiro (^). The former method, while 
it uses the full equation of motion (Equation 1. 7), nevertheless, requires 
some apriori knowledge of the velocity distribution on the axis. In 
addition, the resulting flow field does not necessarily satisfy the 
initial equation in the middle regions. 

Shapiro also discusses a relaxation of finite- difference technique. 
Since this is a purely numerical method that does not lend itself to 
generalized parameter studies, it is not discussed in this document. 

c. Other Methods Developed in Conjunction with the Present Method 

Two approximate methods that originated in the present study are 
of interest because they can be used to check the arbitrarily exact 
solution. 

An investigation of Sauer's method, including the linear coefficient 
Sauer dropped, was undertaken. It was concluded that for a large 
nondimensional radius of curvature ( p), Sauer's method would be 
quite close to an exact solution. However, as p approaches one or 
less than one, Sauer's method does not produce a good approximation. 
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In the vicinity of the nozzle wall, the term 2$y 3?xy» which Sauer 
drops, is not of negligible magnitude. By writing 

= (2 - 14) 


and noting that near the wall 

dy _ 

$ ~ ( 1 + $ x ) M 


dx 


(2.15) 


^ y xy 


9 

dx 


1 _ . 7 

( d y s 

2 

\ 

[ (1 + $x) 2 

' dx 



(2.16) 


where y s = y s (x) denotes the expression for the boundary. 


Thus, by observation of equation 2. 16 it is seen that the condition 
2$y<& X y <<: 1 I s met only f° r small nozzle slopes (large radius of curva- 
tures), particularly on the supersonic side of the throat where (1 $ x ) > 1. 

Accordingly, it should be expected that Sauer's method is invalid for highly 
curved walls, though it should give excellent results for a slowly varying 
nozzle cross-section area. An analogy between the present case and 
external, axisymmetric flow indicates that this is precisely what should be 
expected. In external flow problems, linearized theory is valid only for 
small longitudinal area gradients of the body; therefore, it gives good results 
only when such is the case. 

Method I 

In this method, the procedure used by Sauer is followed with the 
2<$y<I?xy term retained (equation 2. 1). This results in a non-terminating 
series of terms f (x) y 2n , whereas, in Sauer's case all terms are zero 
beyond n = 2. 


The coefficients, assuming fg 
, (7+l) n 2 n- 1 


2 n 


2n 2 


a 


( a n 


= ax, are of the form 
+ b n ax) 


n t 0 

where a n and b n are evaluated from the recursion formulae. 


(2. 17) 


These formulae can be simplified to a recursion involving a n and b n 
rather than f , i. e. , 
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aj = 0 and bj = \ 


while 


a„ = 


b n-l 1 £" 2 a i+l b n _i_i 
4(n-l)* + (i+1) (n-i-1) 


and 


i=n-2 




b i + 1 D n-i-l 
(i+1) (n-i-1) 


n > 2 


(2.18) 


Evaluation of equation 2. 3, using the first five terms of the above 
equations, produces, after differentiation, the following perturbation 
velocity components: 

u = ax + + 1) a 2 y 2 [(7+1) a 2 y 2 l 2 l( 7 +l)a 2 y 2 ] 3 

4 32 144 


and 


11 [(r + 1) a 2 y 2 ] 
6144 


(2. 19) 


_ ( y +l)a‘ 

2 


(y + l) 2 a 3 , 

xy + 16 ; (1 + 2 ax) y 3 


+ (_ L +^_(| + ax) y 5 + 75 g- (T + l)‘<r 7 (i + «*)y 7 


1 


The parameter c can be evaluated by applying the boundary equation 
that the velocity component in the y direction (v) equals zero at y= 1.0 
and x- -e , i. e: 


1-00 


= *v = E < 2i > f 2i (-«> (1.00) 21-1 = 0 

i= 1 


Consequently, 
(7 + 1) 




J=n 


8 


a. + y 2 

1 j=3 J 


(7+ 1) 2j-3 

“i "j 


1+ nii) a .2 + y n Kxtui: 1 g .2(j-D b . 

4 1 L J 1 J 

j=3 


( 2 . 20 ) 


( 2 . 21 ) 


At any point along the nozzle surface, defined by y=y g (x), the 
following should be true*. 


1 _ v 
y s = TRT 

or (2.22) 

V - ( 1 + u) y s ' = 0 

where u is the non-dimensional perturbation velocity in the x direc- 
tion, and y s ' is the slope of the nozzle contour. 

Using 


$y - v and $ x ~ u and substituting equations 2. 3, 2. 17, 2. 18, 
and 2.21 in equation 2.22 results in 


a i ~ 


-Bj + 


}/b* - 4AjCj 

2Ai 


(2.23) 
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t 


where «i is the a value of the i th iteration attempt, and the functions 
Aj , Bi , and Ci are: 



These are solved by starting with an <* value from Sauer's analysis 
and then cycling between the equations until a satisfactory solution is 
achieved. 


In attempting to construct a constant M* line, it should be noted 
from an examination of u and v, as given by equation 2. 19, that 
Sauer's approximation, equation 2. 13, is only good when three terms 
are used to describe the potential function. 
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A better approximation for an M*= 1 curve is 



But, in this analysis (and all further analyses) a curve of an 
arbitrary constant Mach number will be given by 


2* x + *x 2 

+ $y 2 = 

d(2 + d) 

where d (a constant ) is 

defined by M* = 1 + d 

The curve is defined as 

x ~ x(y), where 

LI + 

"v (LI) 2 - 

4(L2) (L0) 


2(L2) 

with 



i= 2n- 1 

j = n-l 


i= n 

L2 = a 1 + ]] 
i= 1 

F' 1 

, 4.,x i+1 2< i+1 > 

(7+1) a 2 

Li 

j=0 

(j+D 

i=n+ 1 

j=i-n 


i= 2n- 1 

j = n-l 


i= n 

LI = Za + Y 

r* 

.i+1 2i+ 1 

(7+1) a 

i= 1 

j=0 

(j+D (i-j) 

i= n+ 1 

j=i-l 


+ b i-j a j+i) y 21 + 

i=n 

y 

, , ,.i 2i+ 1 

(7+1) a 

L 

i= 1 

iZ b i 


(2.25) 


(2.26) 

■s 


. b J ti b i-j v 2i 

(i-j) Y 


>0 


(b j + 1 a i-j + 


2i 


.27a) 
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and 


LO = -d (2+d) + 

i= 2n- 

1 j = n-l 

i= n 

j=i-l 

L l 

y 

i= 1 

LJ 

j=o 

i= n+ 1 

j=i-n 

i= 2n 

j = n 

i=n 

+ 2 

T 

i= 2 

j=i 

i= n+ 1 

j=i-n 

i= n 

+ Y (7+1) 1 or 21 

LJ 

i= 1 



(7+l) i+1 a 2i 

(j + D (H) 


a j+l 



(7+1) 1 a 2l bj bj-j y Zl + 
4j2 (i-j)2 





l (2.27b) 




(NOTE: In using equation 2. 27 when the summation has two sets of 
values, use all the values closest to the summation sign at the same 
time and then use the set of values furthest away.) 


Method II 


Since the interest is in the flow about Mach one, it was theorized 
that a quicker converging series might be obtained if the coordinate 
system was changed so that the axes were the nozzle centerline and the 
Mach one surface. Consequently, this method involves a transformation 
of the coordinate x. 


Let the variable x be replaced by the variable £ , defined so 
that £ = x-g(y), where g(y) denotes the curve representing the critical 
(sonic) line (FIG. 1). The differential equation (including the 2$y4?xy 
term) is transformed from x, y coordinates to the £ , Y frame of 
reference. 
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The g(y) curve is unknown; however, the assumption is made 
that it can be represented as an even polynomial in y, i.e. : 

j=oo 

g(y) = Z a j y 2 ^ (2.28) 

j=° 

where the aj's are, at this point, yet to be resolved. 

Let the x-y origin of the coordinates be located on the axis at the Mach one 
point, in which case a Q = 0. Due to the location of the x-y axis, epsilon is 
still defined as in Sauer's analysis (FIG. 1). The potential function is 
assumed to exist in the form 
i = oo 

* = Z f 2 i U) Y 2i ( 2 - 29 > 

i = 0 


The change of coordinate system produces 


and 




dg 

dy 



-g'f 


> 


(2. 30) 


As in prior analyses, recurrence formulae can be calculated 
by substituting $ and the polynomial form of g(y) into the transformed 
differential equation. This gives f£i (£ ) in terms of f Q (£ ) and the 
aj coefficients. The assumption is made, similar to Sauer's analysis, 
that 


f 0 ' (£)=«£ 

In addition.it is assumed that the coefficients of the Mach one 
polynomial are in the form 

a i = b^ ( 7 +I)* a 2 * * (2.31) 

where bj is not a function of a, x, y, or £ . It becomes convenient 
to represent as 

(7+1) a ^ n ~ 1 t^n + ^ a ^ A (2.32) 

f 2 n 2 n? ” / 
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where A n and Bn are not functions of x, § , or a but may be functions 
of Y or any one of the set of b n 's. 

With the proper manipulation of equation 2. 31 and 2. 32 in the 
simplified equation of motion given by 2. 1, the recurrence formulae 
for A n and Bn are obtained as 

+ 2i 2 bj , when i >. 1 (2. 33) 


and 


A = BiA + 
i 4(i-l)2 


*-i + H- 


ft 


k= i-j-2 


i,i-l 2 


(k+1) 


k= 0 


hk+1 Bj-j-k-ll 

[i-j-k-l)2 J 


(2. 34) 


when i >_ 2 


where 


where 


Ai 0» - 2 


J=1_1 ft ft 

A. = y .gi 
Pl j (i-j) 


, when i >. 2 


(2. 35) 


(2. 36) 


and 


j.i-1 


0 j=i-l 


J x \ 

j*i-l * 
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? 


Before the usual boundary conditions can be applied, the values 
of the bi 's must be found. The restriction to find the necessary extra 
boundary condition arises from the definition of x= g(y) as the curve of 
M* equal one. 

The M* =1.0 curve is necessarily that along which 
2$ x + § x 2 + ®y 2 = 0 

(equation 2.25). The above equations produce a recurrence relation- 
ship for the b^ 1 s, i.e. : 


and 


b 2 = - 


b* = 


1 

“ 4 

1 

32 




J 

j=i-2 

& | 

B J Bi-j 


'*\ 

j 2 (i-j) 2 


(7+1) 

C Ai+1 Ai_i 

l(j+D (i-j) 

k=i-j-2 

- + Z < k+1 > »k + i( 

k=0 \ 

f Ai+i Bi-i_k-l . 
"(j+1) (i-j-k-1) 2 ' 1 ' 

i 


when 


f=i-j-k-2 

j. ^ (l+l) bj+1 Bj Bj-j-k-l -1 


1=0 


j 2 (i-j-k-i - 1) 2 


i > 3 


2. 37 
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The first five terms of the potential equation, when differen- 
tiated, produce velocity perturbation components of 


u 


= a i 


(y+1) 4 a b * (y+1) 5 a 8 8 ■» 

512 y ' 512 y 


and 


_ (Y+1) O’ 2 

* “ 




(Y+1) 2 O' 3 a 

16 y 


> (2. 38) 



while, the first four terms used to describe the Mach one surface are 


s,y) y 2 ( 

(y+4 


1 + or* y 2 + 


— (fl + ») H 


(2. 39) 


(Y+1) 

128 


The boundary conditions that are applied are the same as used 
in Sauer's analysis (equation 2.7). The noticeable difference, however, 
is that while Sauer applies them at x = -t and y=1.0, in this analysis 
the point is t = -|« + g( 1.0)1 and y = 1.0. The resulting equations are 


€ + g(1.0) = - i. 


00 

2 z 

i=0 




(Y+1 ? <* Zi 0 i+ i 


/ (i+D ] 


(2. 40) 
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and 


/> _ , V -L>i + - 4L fcLg 

k ' + * ZLn 


(2.41) 


where 


th 

a is the k iteration attempt to find a 


_ (t+1) 
2 


L, = (p - 1/4) 


I* - 12J1L (p . 1/4, 


= - 1 - z 


i=n i+1 2(i+l) 

(7 + 1) gk-1 


i=2 


(i + 1) 


Bi+1 


2(i+l) 


T ~ P $1+1 


+ Ja±ii [i + (Tr+i) /4] + 


i=l l j=0 J 


l = n r 

Yj (y + l) 1 or 2i /3 i+ ! /(i+1) 
i=0 l k " J 


0,0 " V ( Y+ 1) (p- 1/4) J " (p-1/4) +Vp 2 + 2 P ~ 16 ] 

and n is the number of terms considered in the series. 




2. 42 


(2.43) 


Once the a and c values are known, it is relatively easy to establish 
constant Mach lines. Along constant M*=l + d lines the £ value for any 
y value is given by: 


£ = -X j +^ Xj + X 2 , 


(2.44) 
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where 


i+L k v _1 Ok+l 


. if. V* / ..i+1 2(i +1) 2(i+l)f jii+i . , . v> 

x i =-s-{‘ + .1 <t +1 > * y [ziunr + <t +1 >Z 

j=i-k-l 


k=o 


(k+1) 


( Ai + 1-k 

i + 1-k ‘ 1 
J =o 


(j+l)bj-H /3i-j-k 
(i-j-k)2 


]] 


i +v n j y <7 + D i+2 * 2 < i+1 > 


,-i 


i=0 j= 0 


o + u (i+i-j) 


(2.45) 


and 


= 


d(2+d) 


a 


1 + 


J=i 

_ r, 

i=0 j=0 


OQ Jj” 1 

l 1. 


(7+ 1) 


a 


1 0j +1 

/Si -j+1 

(j+1) 

(i-j+D 


-1 


Methods I and II have a built-in error because they have not taken 
into consideration the $ x 2 and i terms in equation 1. 7, but they 
have included similar order of magnitude terms in, the series expansion 
of $ . 


Neither methods I nor II have yet been programmed or substituted 
to replace Sauer’s analysis. But this will be done later, and further 
comparison and analysis will be attempted. 

3. PRESENT METHOD 

In view of the requirements for an analysis giving exact results, it 
was considered advantageous to proceed directly to a solution of the 
exact equation of motion. For convenience, that equation is repeated 
here. 

[l' 2 + *x)®x - r ®y 2 ] ®xx + [( 7 TT) • 1 ®x< 2 *»x>-®y 2 j® yy - 

" TyTT) (1+ ®*> ®y ®xy + (l7TTT ‘ r [®* (24 ? x) ^y 2 ]^ = o(3.i) 

The potential function is assumed to exist in the form 
i=o° j =oo 

® = Y. Y A(i, 2j)x i y 2 J ( 3 - 2) 

i=° j=0 
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where the even exponent of the variable y is caused by the symmetry of 
the flow field. The velocity distribution along the x axis would be 

y i A(i,0)x i 1 
i=l 

and not just the first approximation used by Sauer (o-x) 

Upon differentiation and substitution of equation 3. 2 into 3. 1, two 
algebraic relationships result. These relations, which are the necessary 
and sufficient conditions for the assumed potential function to satisfy the 
equation of motion, are the recurrence formulae defining the A(i, 2j) 
coefficients in terms of the set of A(i, o) coefficients. The recurrence 
formulae result from the existence condition that the coefficient of each 
x and y term in the series must vanish separately in order to satisfy 
the equation of motion 3.1, for an arbitrary x and y value. 


These recurrence formulae are given below. Their formidable 
appearance is the direct result of the extreme nonlinearity of the equa- 
tion of motion. With the definition of 



( 3 . 3 ) 


2(k+2)(k+l)(i+l-k)A(k+2, 0)A(i+l-k, 0) + 


+ 8 T (i+l-k)A(k, 2)A(i +l-k, 0) + 
a=i-k- 1 

+ E °a o [(k+2) (k+1) (a+1) (i+l-a-k) A (k+2, 0) A(a+ 1 , 0) A(i+l-a-k, 0) + 

a=0 

+ 4 I (a+1) (i+l-a-k) A(k, 2) A(a+1, 0) A(i+l-a-k, 0)] } 


( 3 . 4 ) 
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and 


r . k=i- 1 

A(i,2(j+2» = 8 - jp^ r P'lA 8(i+1 ' k)( j +2, * r A ( k ' 2(j+2 ^ A(i+1 " k * 0)+ 

+ 2(k+2) (k+1) (i+l-k) A (k+2, 2(j+l))A(i+l-k, 0)+ 
a=i-k-l 

+ Yj [(a+1) (i+l-a-k) (k+2) (k+1) A (k+2, 2(j+l))A(a+l, 0)A(i+l-a-k, 0) + 
a=0 


+ (a+1) (i+l-a-k)4(j+2) 2 T A(k, 2(j+2))A(a+l, 0)A(i+l-a-k, 0)] J + 
k=i n=j 

+ Z Z I 2(i+l-k)(k+2)(k+l)A(k+2,2n)A(i+l-k,2(j+l-n) + 

k=0 n=0 1 


+ (i+1 -k)8(n+l) 2 T A(k, 2(n+l))A(i+l-k, 2(j+l-n)) + 


+ (y-p) (n + l)(k+l)(j +l-n)A(k+l, 2(n+l))A(i-k, 2(j+l-n)) • 
a=i-k 

+ Y f (a+l)(i+l-a^kH(n+l) 2 r A(k, 2(n+l))A(a+l , 2(j+l-n))A(i+l-a-k, 0) + 
a=0 l 

+ (a+l)(L+l-a-k)(k+2)(k+l)A(k+2, 2n)A(a+l, 2(j+l-n))A( i+l-a-k, 0) + 
b=n-n 

+ V {(a+l)(i+l-a-k)4(n+l) 2 T A(k, 2(n+l))A(a+l, 2b)A(i+l-a-k, 2(j+l-b-n)) + 
b=0 

+ (a+l)(i+l-a-k) (k+2)(k+l)A(k+2, 2n)A(a+l, 2t)A.(i+l-a-k, 2(j+l-b-n)) + 

+ (j+l-n-b)(k+2)(k+l)(b+l)4 T A(k+2, 2n)A(a. 2(b+l))A(i-a-k, 2(j+l-b-n) + 


+ (j+l-n-b)(n+l)(k+l)(a+l)(y^) A(k+1, 2(n+l))A(a+l, 2b)A(i-k-a, 2(j+l-n-b)) + 
+ (j+l-n-b)(n+l)(b+l) (16) ^n+1 - yjy) x 

(k, 2(n+l))A(a, 2(b+l))A(i-a-k, 2(j+l-b-n))}] J j 


x A 
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The original written outline of a method to handle the programming 
of the two recurrence relations (equations 3.4 and 3. 5) had the Ai,2j's 
resolved into polynomials in the unknowns o' 8, All advice obtained 
suggested that the Ai, 2j's be evaluated as pure numbers instead of 
functional representations of y and the unknowns, this suggestion was 
followed, but one third of the computer running time will be saved if 
the program is revised to its original form. 


For convenience in future discussions, let the unknowns be repre- 
sented by D^'s , where 


D i = A i+ 1 , 0 


(3.6) 


when 


i > 1 


The equations are now being re-programmed to appear in the form 
of polynomials of the type 


i=k 

A i. 2j = Z e i D ! 

i=l 


1, i d 2,i d n, i d 0, i d p , i 


D 


-- D 


n 


, (3.7) 


where 

ej[ is a pure number, d n ,i is the power to which the D n unknown is 
raised in the i^ term 


and 


0 = 7+1 

As equation 3. 4 is already programmed in the above form, some of the 
Ai, 2 s are presented in Table I. 

Here, as in Sauer's analysis, the origin of coordinates is chosen 
on the axis a distance e downstream of the minimum section. 

At this point in the analysis the unknown functions have been reduced 
to a system of unknown coefficients of the type Di and also the unknown « . 
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The general inviscid flow boundary condition is that the velocity 
at the boundary surface must be tangent to that boundary. The boundary 
surface y g is represented as a power series in x 
i= oo 

y s = Z B i xi < 3 * 8 > 

i=0 

While, the inviscid boundary condition can be written as 
(1+u) - v = 0 


where both u and v are evaluated along the boundary surface, i. e: 


u = u(y s , x) = u(£B[X 1 ,x) = u s (x only) 


and 


(3. 10) 


v = v(y g , x) = v(Z 


Bjx 1 , x) 


vs(x only) 


Substituting equations 3. 10 and 3. 8 into 3. 9 and collecting the 
resulting terms as a power series in x, results in the boundary condi- 
tion equation in the form of 
i=oo . . 

2 G(i,w) x =0 (3.11) 

i=l 


where 


j=oo q=i 

' 

a=i - q 

G(i,oo) = iB^ + Yj Yj ' 


Yj (a+l)q B(a+1) C(i-q-a, 2j) 

ii 

o 

+3 

II 

O 


a = 0 


-2 <r q . (j)C(i-l-q, 2j-l)| A(q, 2j) (3.12) 

with C (i, n) given by 


C(i, 0) = j+1 for i = 0 | 

l 0 for i^= 0 j 


(3. 13) 
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and 


C(i, 1) = Bi 


1 


C(i, n) 


b=i 


Z 


B, C. , . 

b i -b, n- 1 


with 


n > 2 


(3. 14) 


In equation 3. 11, the G(i,oo) designation simply denotes the total 
collection of terms constituting the coefficient of x 1- ^. G(i,oo) con- 
tains an infinite number of columns of terms nonlinear in Ap , ^ and 
hence the oo notation. 


The only possible solution of equation (3. 11) that would be valid 
for any arbitrary value of x, would be that 


G(i, oo ) = 0 


(3. 15) 


for all values of i (i=l, 2, 3 oo). 

Equation (3. 15) thus denotes an infinite set of nonlinear simultaneous 

equations in the unknowns (i=l, 2, 3 oo), each one containing an 

infinite number of terms. 

In computing numerical solutions the infinite set of equations 
need not be considered, but only those equations that will insure con- 
vergence to the required accuracy in the domain of interest. Further- 
more, only the number of terms in each equation need be considered that 
are necessary to produce the required accuracy. 



(3. 16) 
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Remembering that the first column constitutes the unknowns with 
all the other coefficients expressible in terms of the first column, 
(through application of equations 3. 4 and 3. 5) a maximum number of 
rows, M+l, can be considered, together with a maximum number of 
columns, p+1. All rows and columns beyond M+l and p+1 can be 
considered negligible. 


The expression G(i,oo) = 0 (equation 3. 15) contains all terms 
occurring in the first i+1 rows of 3. 16 and an infinite number of its 
columns. For the M+l, p+1 system the expression G(i,oo) is changed 
in notation to read G(m, p), thus, 


G(m, p) = 0 

where 

m = 1, 2, 3 M 


(3. 17) 


The expression for G(m, p) is the same as equation .3. 12 with i re- 
placed by m and oo replaced by p. This system (equation 3. 17) has 
M equations and contains the terms of the first M+l rows and p+1 
columns of the Aj^ 2j array, 3. 16. 

In using the recurrence equations, 3. 4 and 3. 5, to evaluate the 
Ai, 2j coefficients (i<M, j<.P) considered in the boundary condition 
equations (G m p's), certain A^ Q terms for which k >M+1 appear. This 
is due to diagonal propagation characteristics of the unknowns through- 
out the array of 3. 16. These unknowns must be considered zero in the 
array of coefficients i< M, j< p. This consideration is justifiable on 
the grounds of the expected small magnitudes of these unknowns. 

A judicious initial choice for M and p values must necessarily 
result from experience. This choice shall be discussed in Section VI b. 

The system is now reduced to one containing M equations and M 
number of unknown A^ Q 's. The quantity e, which locates the origin 
of coordinates, is an additional unknown, giving a total of M+l unknowns. 
The M+lst constraint, necessary to evaluate e, stems from mass flow 
considerations and will be discussed later. 

At this point in the analysis the major task remaining is the 

evaluation of the unknowns £ & D m . Calculation of these 

quantities allows the calculation of the velocities, u and v, and hence 
the calculation of all other flow properties. The terms are to be 
calculated for various prescribed e values. An iterative scheme 
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necessary to compute the D{ terms is discussed below. 


The procedure used in the present computer program is as follows: 

The i^ equation will be used to solve for the i^ unknown by arranging 
G(i,p) in the form (note: there is no difference between G^p or G(i,p)) 

Gi, p = ai ; 0(Di,D 2 --- Di_ 1; D i+1 -- > D M )+ai > i(D 1> D 2t --D i _ 1( D i+1 --,D M )D i + 

+a i, 2i D l, D 2, "~ D i-l, D i+l > D M) D i 2+ 

+ a i, n( D l> d 2, "“ D i-1, D i+1"'> D M) D^+etc. (3.18) 

where Dj represents the unknown of G(i,p) = 0. Equation 3.18 is 
simply an expression for G(i, p) written as a polynominal in D(i) where 
the coefficients are functions of Dl, D 2 °M exclusive of D]_. 

The present program used three values of Dj (holding all others 
constant) to fit a cubic formula to vs G^ p. Using the cubic formula, 
a G{ p = 0 point is predicted for a certain D| value. After that Dj value 
is tried, a more accurate cubic is fitted; this continued until an actual 
Gi, p = 0 point is found or 50 tries are made. 

The process is repeated for all values of i in equation 3.18 

(i = 1, 2 M). A suitable root for each G^ p is used as new input in 

the Gi-fi p equation to generate a second set of coefficients and so on - 
until the final equation (Gj^p). After the process is performed for the 
equation, all of the G(i,p)'s are evaluated with the on hand values of 
the D's. If the absolute values of each and every one of the boundary 
condition equations is less than or equal to a specified tolerance ( a), 
then the on-hand D values are considered a good solution. 


Let a cycle be defined as a set of equation solutions - from i=l to 
i=M of G(i,p) - and let k be the number of cycles performed. After the 
initial cycle (k=0), G m) p is apriori non-zeroed. To explain, let G m, p 
be equal to minus a number (n) times the value of the G m> p from the 
k-1 iteration cycle, i. e. , 


m, p 


= -n G. 


m, p 


k-1 


The problem is to determine n. 


(3.19) 
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On the last equation (M^), it was necessary to leave n equal to 
zero, as no perturbation occurred to this equation after a Dj^ is found. 
In all the other equations, on the k=l and k-Z cycle it was determined 
best to let n=. 75. After the k=2nd cycle, n could become a function 
of the equation as well as the cycle, i. e. 


Gm.pj _ u G r 


m, k ^m, p| 


(3. 20) 


k-1 


where 


n m, k =n m, k-1 + q (3. 21) 

On the mth equation and the k^h cycle a q is calculated from 


G 


m,p 


q = 


k-1 


(3. 22) 


m.p 


k-2 


With the restrictions on equation 3. 21 as follows 
If 3_<c_|q|< 50, then n m> k =n m> k _j ± 0. 3 0 
depending on the sign of q 


Or if 


q^ 50, then n 




m, k l m, k-1 


(3. 23) 


(3. 24) 


While /3 is an input parameter to the program, it was found best to 
let it equal 0. 1 or zero. 


The a priori non-zeroing procedure is repeated until every 



a 


m = 1 , 2 M 


(3. 25) 


It is obvious from the method just described that the closer the 
initial estimated D values are to the actual solution D values, the 
quicker the computer will obtain the solution. Consequently, it 
becomes advantageous to prepare a library of solutions that can 
be consulted before a new problem is started. The library obtained 
appears in the Appendix. 
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The remaining variable is e . Neither a one dimensional analysis 
nor Sauer's analysis allows a variation in mass flow in the nozzle. As 
it is physically possible to have a slightly changing nondimensionalized 
mass flow through a rocket nozzle, it should be possible, in theoretical 
calculations, to have varying mass flow conditions without violating the 
theoretical restrictions on the problem. Consequently, it should be 
possible to allow a variation of e , so as to allow a variation in the 
nozzle mass flow. 

The variation in € will be physically restricted to a range of values 
between choking conditions and no supersonic flow. Violation of the 
physical restrictions in e might show up in the numerical analysis. 

There are two different ways in which too high, or too low an 
value will appear. 

Physically, coo high an e value can show up as a supersonic flow 
stream with a subsonic pocket ({FIG 2) . What seems to happen is 
that as e increases the mass flow also increases to the point where it 
is impossible to have a completely subsonic cross-section. Consequently, 
the Mach one line falls back to the axis instead of proceeding to the nozzle 
wall. With too low an e value it is expected that the Mach one line will 
intersect the nozzle downstream of the minimum section, A downstream 
intersection would probably violate the entropy law. 

One of the mathematic assumptions made quite early in the analysis 
was that 


D i+1 


« Di 


(3. 26) 


for any i value. Consequently, when the absolute value of D 2 approaches 
the value of D]_, it is necessary to question the mathematical validity of 
the resulting solution (this occurs when € is too large or too small.) 
Applying a rule of thumb: A mathematically valid solution shall be 
considered one where 



< . 20 


(3.27) 


The above statement, necessarily, restricts the range of € . In 

the problems observed so far, the lower limit of e has always been 
Stipulated by the mathematical limitations, while the upper limit has been 
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NOIllSOd A 



FIGURE 2 SUPERSONIC FLOW WITH SUBSONIC POCKET 


stipulated either by the mathematic limitation or by the physical 
limitation. 

Since e is defined implicitly in terms of mass flow, it is necessary 
to solve the problem for various values of e and then calculate the 
mass flow for each of these solutions. The e values are then plotted 
versus mass flow. If the mass flow is prescribed, the corresponding 
e is read from the plot, and a final solution can be obtained on that basis. 
A tentative plot of e versus non dimensionzlied mass flow appears in 
the Appendix as FIG A-9. 

4. ANALYSIS AND USE OF THE DATA 
FROM THE PRESENT METHOD 

Once a solution is found, the problem has just begun. From the 
definition of the potential function, equation 3. Z, the perturbation 


( i +1 )A(i+l , Zj) x 1 y 2 J (4. 1) 


Z(j + l)A(i, Z(j+1) x 1 y 2 i + 1 (4. Z) 


All of the A^ ; zj ' s are evaluated as pure numbers for the particular 
flow field for which a solution exists. Consequently, at any point in 
the flow field the Mach number, velocity components, and non-dimension- 
alized pressure and density can be evaluated. 

If M* is defined in terms of a parameter dj as 

M* = l+d x (4.3) 

then the equation 

Zu + u 2 + v 2 - d,(Z + d,) =0 (4.4) 


velocity components u and v are: 

i=M j=p 
u = u(x, y) = £ £ 

i=0 j = 0 


and 


i=M j=p 

v=v(x,y)= 2 

i=0 j = 0 
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describes a constant M* curve (where M* is given by equation 4. 3). 
Along the curve the velocity components and the flow angle (#) are 
evaluated, the velocity components by equations 4. 1 and 4. 2 and the 
flow angle by 

0 =tan ’ l T^ < 4 - 5 > 


The initial equation to be satisfied was equation 3. 1. An idea can 
be obtained as to whether a point in the flow field is within the region 
of convergence of the solution or not by evaluating equation 3. 1 to see 
how close it comes to zero. Define that number as zero. 


At the intersection of the constant curve and the nozzle, a 
check is made as to how close the solution satisfies the imposed 
boundary condition. Define the error of the difference between the 
stream slope and the nozzle wall slope as er, where 


er 



v 

1+u * 


y S ' 


(4.6) 


The constant M* line is used as input in the characteristics 
nozzle design program. Which M* line is chosen depends on. the 
error printout values (zero and er). 

Another type of output available is a straight line and is defined as 
starting on the axis at a given position and extending to x=0 and the 
nozzle wall. 


At each increment in position along such a line the flow angle, Mach 
number and M* value are obtained. 


Along the nozzle surface the M*, y s ', tan 6, Zero, and er values 
are calculated. 


With all the available printouts, it is not too difficult a problem 
to analyze the resulting flow field and stipulate the region of applica- 
bility of the solution. 
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5. INPUT PARAMETERS 


The present method assumes that the minimum section of the 
nozzle's surface can be expressed as a power series in x, i.e,.’ 


i = oo 

y s = 2 B i xi (5-D 

i = o 

(a repeat of equation 3. 8). 


Before a problem can be started for a given nozzle contour, a power 
series must be fitted to the nozzle surface. If the minimum section can 
be represented as an arc of a circle, one can use the equation 

(x+O 2 + Jy - (l+p)j 2 = p 2 (5.2) 

where p and e are defined as shown on FIG 1. 


The above equation can be written as 


y = f(x) = 1+p 




A MacLaurin series can then be applied in the form 


(5. 3) 


y = f(0) + f' (0) x + 


f"(o) 

”21 


X 


+"( 0)1 

PH 


etc, 


(5.4) 


where f ( 0) , f'(0), and so on, can be evaluated from equation 5. 3. 

After differentiating equation 5.3 the necessary number of times, 
the values are found as 


B 0 = 1 + P (1 - g) 


(5. 5) 
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One of the purposes of this analysis is to be able to handle a non- 
constant radius of curvature solution. A subprogram to evaluate any 
nozzle minimum section in terms of a power series is in process. 

DISCUSSION 

1. MATHEMATICAL PROBLEM AREAS OF THE PROGRAM 

Some mathematical areas to be clarified from the preliminary 
report (Ref 5) are: 

a. Non-Unique Solutions 


The form of the boundary condition equations (equations 3. 12 
and 3. 17) is such that more than one solution for is possible in the 
iteration scheme. While a certain number of them will be complex 
conjugates, more than one might be real. There are two possible 
physical explanations for more than one real solution. If a particular 
solution of a particular boundary equation would not iterate success- 
fully, one might be investigating an unstable solution that could not 
exist physically for more than an instant of time. Needless to say, in 
the transonic flow region, unstable solutions have become infamous. 
Even if more than one complete problem solution set is found, it is 
suspected that all except one would represent a diverging series. A 
diverging series could quite possibly represent another unstable 
solution or possibly be the result of extraneous roots. A diverging 
solution, mathematically, cannot be considered a valid solution. 

In most cases, the starting data represents a solution similar 
to the one desired. Consequently, it is felt that the iteration scheme 
solves for the closest one which is the right one. Also, if the program 
tries to come up with an extraneous solution, the iteration scheme 
might oscillate until the program found itself back at the valid case. 

While the constant M* lines are being constructed, the multi - 
value problem arises again. For any given y, there is more than one 
x value that lies on the constant M* curve because of the nature of 
u and v, and conversely. 

Experience with the program has shown that the transonic range 
is repeating itself, just as a sine or cosine wave does every 360°. The 
region of interest, however, is the vicinity of the throat of the nozzle. 
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Consequently, the program is looking for the particular curve 
that starts at the x axis, close to the origin, and propagates upward, 
toward y greater than one. 

The program for calculating the M* lines is written so it 
will reject any point it finds that is further away from the previous 
point than four times the step size but will find the right point instead. 

b. Region of Applicability 

Let the region of applicability be defined as that portion of the 
flow field where the boundary condition and initial equation are satisfied 
to a desired degree of accuracy. The number of unknowns that need 
be considered is directly a function of the region of interest. If one 
was interested only in the region about the y axis (x quite small), then a 
small number of unknowns need be considered. On the other hand, the 
larger the highest value of x becomes, the more unknowns that are 
required to satisfy a given tolerance. 

It is usual in transonic analyses to always include, inside 
the region of applicability, the Mach one line. But, in this analysis, to 
include the Mach one line necessitates use of quite a few more unknowns 
than would otherwise be needed. Consequently, with the proper restric- 
tions, the analysis can be limited to four unknowns. 

The constant M* line that would become the output of the pro- 
gram must be chosen prudently by analyzing the various error functions. 
Fortunately.it turns out that as the radius of curvature becomes 
smaller the preferred M* value becomes higher, producing better 
accuracy in the next analysis (characteristics). 

In the case of a radius of curvature of 4. 0, a plot of various 
M* lines for a M, p system of 4, 7 and 7, 9 (FIG 3) shows that 
the M* = l. 0 lines vary considerably. But, the higher valued M* lines 
for the 4 and 7 unknown systems are similar (M*=l. 10 through 1. 15), 
with the lines closest to x=0, y=y g (M* = l. 13) practically identical. 
Consequently, if the M* = l. 13 line was used as the output data, only a 
4 unknown system would be required. Along the M* = 1. 13 line for 
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FOR A MINIMUM SECTION RADIUS OF CURVATURE OF FOUR 









the 4 unknown case, the maximum zero value is . 001 with er = + . 00005. 
Therefore, M*=l. 13 is a good choice of data to be used in the charac- 
teristic program. 

The choice of the number of columns to use is more difficult 
than the choice of the number of unknowns. The main restriction to 
be satisfied is that the last column should have no effect on the region 
of interest. The problem complicates itself if too many columns are 
included because too many possible solutions are allowed. At a certain 
point the addition of another column introduces a solution that is slightly 
different than the solution required, leaving the computer with no way 
of choosing which solution is better. Experience has shown that 9 
columns is more than sufficient but yet not too numerous to cause 
extraneous solutions. 


The array of ' s for a p = 2.00 and e = .165 (Table II) 

shows that the effect of the ninth column for an x value less than 
± . 2 would be small compared to the other columns. 


that 


One of the initial assumptions made about the Aj £j series is 


lim 
i -*• oo 


and 


lim 
j— oo 


A i, 2j X " 


j = l, 2- 


-oo 


L i-2j 


< a 


( 6 . 1 ) 


i=l, 2- 


- 00 , 


in all cases checked the above has been true. 

One of the minor problems in handling the program was the 
discovery of unbalanced arrays of A^ An unbalanced array 

can be defined as that array where the overemphasis on one of the 
variables produces a solution that has no physical significance. With 
the case of M=4 and p=3 the x variable is overemphasized. This occurs 
when the number of columns (p+1) is smaller than the number of 
rows (M+l). There exists the implicit assumption that the x variable 
exerts a greater effect than the y variable. However in the actual 
physical system, the converse is true; the boundary condition is applied 
where y>l. 0 and x<l. 0. Consequently, on physical grounds alone, one 
is justified in not trusting any solution where M>p. 
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2. UTILIZATION OF THE PRESENT PROGRAM 


This analysis is just one of three parts of a procedure to design 
a given rocket nozzle. Initially, the propellants, chamber pressure, 
fuel to oxidizer ratio, minimum section contour radius, nozzle area ratio, 
and the rocket nozzle length are specified. A chemical equilibrium program 
is used to evaluate the chemical reactions and find the pressure and ratio 
of specific heats (7) at various Mach numbers. The present analysis 
utilizes the equilibrium 7 value to produce input data for a method of 
characteristics equilibrium program that then designs the diverging 
part of the nozzle for various exit pressures. 

Using equilibrium data as an input parameter for a frozen flow 
program (Transonic Analysis) and then using these results in an equilib- 
rium program (Characteristics Program) seems contradictory. 


Although the total flow field is to be treated as an equilibrium 
problem, one part of that field (under certain conditions) can be treated 
as frozen flow if 


where 6 is the significant length parameter of the frozen flow region, 

(FIG 4 ), Y is the lowest velocity through that region, and R t is 
the fastest reaction rate of the chemical processes that might occur in 
the region. As Y and R t cannot be varied in the programs, the only 
choice is to reduce, as much as possible the size of 6. The region defined 
by 6 is reduced by using, in the transonic analysis, the 7 value from the 
chemical equilibrium program which corresponds, as closely as possible, 
to the Mach number of the most accurate constant Mach number line. 


3. COMPARISON WITH OTHER ANALYSIS 


It was predicted previously that this analysis and Sauer's would 
tend to, agree more with one another as the radius of curvature increases. 
For the case of a radius of curvature of five, there is already a noticeable 
difference as can be seen in FIG 5. 


There is a band of Mach numbers where a transonic analysis 
and the method of characteristics are both applicable. If the transonic 
data agreed with the characteristic data in this band one could draw the 
conclusion that the transonic method used was a valid one. 


40 


to 



41 


FIGURE 4 FLOW FIELD OF THE THREE PROGRAMS 




BY SAUER’S METHOD AND THE PRESENT ANALYSIS 




A comparison was made utilizing a nozzle of radius of curvature 
of two. First, two different M* lines were calculated by Sauer's method 
as inputs for the characteristic program. The characteristic program 
proceeded, in each case, to calculate the Mach number distribution 
along the nozzle wall (FIG 6). The two different curves are symbolized 
by a0 and a A . From Sauer's equations, the hypothesized Mach 
number distribution was calculated - it is represented by on 
FIG 6. The characteristic curves and the predicted curve do not 
agree too well with each other. 

The procedure was repeated with the present transonic analysis 
calculating the input M* line for the characteristic program. The 
resulting Mach number distribution is symbolized by © on FIG 6. It 
can be seen on the figure that the Mach number distribution along the 
wall calculated using the present method (©) lies quite close to the 
one calculated by the method of characteristics (©). Consequently, the 
present transonic analysis seems to be accurate. 


RECOMMENDATIONS 

Before the conclusion of this study, the following has yet 
to be completed: 

1. SIMPLIFIED ANALYSES 

The two simplified analyses (Methods I and II) must be tried 
for various cases. Most important is deciding how many terms are 
to be used in their potential functions infinite series. After numerical 
cases are tried, both methods will be used in place of Sauer's method 
in the Allison program (6). Of course, a comparison between Sauer's 
analysis. Method I, Method II, and the present method will be made. 

2. LITERATURE SEARCH 

The literature search, now in progress, will be continued 
until all available report, books, and periodicals have been reviewed 
for useful methods and experimental data. 
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MACH NUMBER ALONG THE WALL OF THE NOZZLE 
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FIGURE 6 COMPARISON OF PREDICTED MACH NUMBER 
DISTRIBUTION WITH ACTUAL DISTRIBUTION 
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3. PRESENT METHOD 


The Ai t 2j sub -<rou tine must be re -programmed as mentioned in 
Section HI. Afterwards, it is expected that other changes will become 
apparent that can further shorten the running time of the program. 

After these changes are made, the p value will be pushed as 
low as the program will go. The resulting listing of all the unknowns, 
for particular p and e values, will become part of the program for 
use with a routine for choosing initial values of the unknowns to start 
the program. Most of the solutions obtained are based on a 7 of 1. 20; 
to efficiently handle other 7 values further work will have to be done. 
Increments of the unknowns, (not necessarily small) versus increments 
in 7 values, must be obtained for various p and € values. 

The nozzle parameter sub-routine will soon be revised so it 
can handle non-constant radius of curvature cases. After no more 
program changes are necessary the programmer will streamline 
and polish it up, so that anyone can utilize it with a minimum of 
instructions. Such instructions will appear in the final report. 

4. FINAL REPORT 

The final report will contain all of the material presented herein, 
including revisions and changes. Examples and trends discovered 
in the course of the evaluation of the analysis will be included. 
Derivations of the equations used, as well as the operational methods, 
will appear as appendixes. 

The report will contain a write-up of all the programs, with 
a section written by the programmer presenting symbolic write-ups, 
instructions for use, and flow charts of all the sub-programs in use. 
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APPENDIX 


LIBRARY OF SOLUTIONS 


The solutions achieved to date appear in two forms, the unknowns 
versus the radius of curvature of the minimum section, p, and the 
unknowns versus the distance, e, of the Mach one line from the minimum 
section. 

The last curve (FIG A-9) represents the various non dimensionalized 
mass flows (for .constant p values) versus e values. The relationship 
between the mass flow (m) and the non dimensionalized mass flow (m) 
is 

0 ✓ — ' 

m = y2 p a* m, 
o o 

where p Q is the isentropic stagnation density. 
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Coeff, 



TABLE I 

2 CONSIDERING ONLY 4 UNKNOWNS 



A 2, 2 
471 ) 
1 . 0 
2 . 0 


1 1 

3 0 

3 0 


0 0 
0 0 
0 0 


1 0 
1 0 
2 1 


A 3, 2 
8 . 0 
6 . 0 
4. 5 
12. 0 
4. 0 
4. 0 


1 

2 

0 

2 

4 

4 


0 

1 

2 

1 

0 

0 


1 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 


1 0 
1 0 


1 

2 

2 

3 


0 

1 

1 

2 


A 4, 2 
12. 5 
10 . 0 
11.25 
15. 0 
20 . 0 
30. 0 
22 . 5 
30. 0 
2 . 0 
12.0 
8 . 0 


1 

2 

1 

0 

2 

3 

1 

3 

5 

5 

5 


0 

0 

2 

1 

0 

1 

2 

1 

0 

0 

0 


0 

1 

0 

1 

1 

0 

0 

0 

0 

0 

0 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 

1 

1 

1 

2 

2 

2 

3 

2 

3 

4 


0 

0 

0 

0 

1 

1 

1 

2 

1 

2 

3 


a 5 , 2 
15. 0 
36. 0 
22. 5 
30. 0 
48. 0 


2 

1 

0 

2 

3 


0 

1 

1 

0 

0 


0 

1 

0 

0 

1 


1 

0 

1 

1 

0 


1 0 
1 0 
1 0 
2 1 
2 1 



Coeff. 
(Cont. ) 

A 5, 2 
6. 75 
81. 0 
12. 0 
72. 0 
48. 0 
18. 0 
108. 0 
13. 5 
81. 0 
72. 0 
12 . 0 
32. 0 
16 . 0 

a 6 , 2 

28. 0 
70. 0 
31. 5 
252. 0 
35. 0 
105. 0 
70. 0 
28. 0 
168 . 0 
94. 5 
63. 0 
378. 0 
56. 0 
63. 0 
252. 0 
112 . 0 
126 . 0 


TABLE I (CONTINUED) 

A i( 2 CONSIDERING ONLY 4 UNKNOWNS 


Power of: 

D^ D£ D3 D^ y r 


0 

2 

0 

1 

3 

4 
4 
0 
2 
4 
6 
6 
6 


3 

2 

0 

1 

0 

1 

1 

3 

2 

1 

0 

0 

0 


0 

0 

2 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 

2 

1 

2 

3 

2 

3 

2 

3 

4 

3 

4 

5 


0 

1 

0 

1 

2 

1 

2 

1 

2 

3 

2 

3 

4 


1 

1 

3 

0 

2 

0 

1 

3 

4 
4 
1 
3 

3 
1 
0 
2 

4 

5 


1 

0 

0 

2 

1 

0 

1 

0 

0 

0 

3 

2 

2 

0 

2 

1 

0 

1 


0 

2 

0 

1 

1 

1 

0 

0 

1 

1 

0 

0 

0 

2 

1 

1 

1 

0 


1 

0 

1 

0 

0 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


1 

1 

2 

1 

2 

1 

2 

3 

2 

3 

2 

2 

3 

2 

2 

3 

4 
3 


0 

0 

1 

0 

1 

0 

1 

2 

1 

2 

1 

1 

2 

1 

1 

2 

3 

2 
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TABLE I (CONTINUED) 

2 CONSIDERING ONLY 4 UNKNOWNS 


Coeff. 
(Cont. ) 

A 6, 2 
336. 0 
94. 5 
252. 0 
168. 0 
4. 0 
48. 0 
80. 0 
32. 0 


Power of: 


Di 


D: 


D- 


D/ 


7 + 1 


5 

1 

3 

5 

7 

7 

7 

7 


1 0 0 

3 0 0 

2 0 0 

1 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 


4 

3 

4 

5 

3 

4 

5 

6 


3 

2 

3 

4 
2 

3 

4 

5 


A 7, 2 
80.0 

1 

0 

1 

1 

45.0 

0 

2 

0 

1 

36. 0 

2 

1 

0 

1 

48.0 

0 

1 

2 

0 

192.0 

2 

0 

2 

0 

40. 0 

4 

0 

0 

1 

240. 0 

4 

0 

0 

1 

432. 0 

1 

2 

1 

0 

192. 0 

3 

1 

1 

0 

1152. 0 

3 

1 

1 

0 

25. 0 

.0 

0 

0 

2 

160. 0 

1 

0 

1 

1 

90.0 

0 

2 

0 

1 

360. 0 

2 

1 

0 

1 

160. 0 

4 

0 

0 

1 

192. 0 

5 

0 

1 

0 

512. 0 

5 

0 

1 

0 

108. 0 

2 

3 

0 

0 

40. 5 

0 

4 

0 

0 

648. 0 

2 

3 

0 

0 

540. 0 

4 

2 

0 

0 

1440.0 

4 

2 

0 

0 


1 

1 

2 

1 

2 

2 

3 

2 

2 

3 

1 

2 

2 

3 

4 

3 

4 
2 
2 
3 

3 

4 


0 

0 

1 

0 

1 

1 

2 

1 

1 

2 

0 

1 

1 

2 

3 

2 

3 

1 

1 

2 

2 

3 




TABLE I (CONCLUDED) 

A i( 2 CONSIDERING ONLY 4 UNKNOWNS 


Coeff. 

( Cone. ) 


a 7, 2 

96 . 0 
192 . 0 
432. 0 
768. 0 
256. 0 
48. 0 
576. 0 
960 . 0 
40. 5 
432. 0 
720. 0 
384. 0 
32. 0 
160 . 0 
192. 0 
64. 0 


Power of: 

Dj d^ D3 D4 7+1 r 


0 

2 

1 

3 

5 

6 
6 
6 
0 
2 

4 
6 
8 
8 
8 
8 


1 

0 

2 

1 

0 

1 

1 

1 

4 

3 

2 

1 

0 

0 

0 

0 


2 

2 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


2 

3 

3 

4 

5 

3 

4 

5 

3 

4 

5 

6 

4 

5 

6 

7 


1 

2 

2 

3 

4 
2 

3 

4 
2 

3 

4 

5 

3 

4 

5 

6 
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TABLE II 
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FIRST UNKNOWN 



P RADIUS OF CURVATURE OF MINIMUM SECTION 


FIGURE A1 
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FIRST UNKNOWN VERSUS RADIUS OF CURVATURE OF MINIMUM 
SECTION FOR CONSTANT VALUES OF THE DISTANCE BETWEEN 
MINIMUM SECTION AND MACH ONE LINE 
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FIGURE A3 THIRD UNKNOWN VERSUS RADIUS OF CURVATURE OF MINIMUM SECTION FOR 
CONSTANT VALUES OF THE DISTANCE BETWEEN MACH ONE LINE AND MINIMUM SECTION 













.10 .11 .12 .13 .14 .15 .16 .17 .18 .19 

€ ; DISTANCE FROM MINIMUM SECTION TO MACH ONE LINE 


FIGURE A5 FIRST UNKNOWN VERSUS DISTANCE FROM MINIMUM SECTION TO MACH 
ONE LINE FOR CONSTANT VALUES OF RADIUS OF CURVATURE OF MINIMUM SECTION 
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_ 04 i TrrnrTTT f TTTTTn rnTrrrnTTF 

£j DISTANCE FROM MINIMUM SECTION TO MACH ONE LINE 


FIGURE A6 SECOND UNKNOWN VERSUS DISTANCE FROM MINIMUM SECTION TO MACH 
ONE LINE FOR CONSTANT VALUES OF RADIUS OF CURVATURE OF MINIMUM SECTION 
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FIGURE A7 THIRD UNKNOWN VERSUS DISTANCE FROM MINIMUM SECTION TO MACH 
ONE LINE FOR CONSTANT VALUES OF RADIUS OF CURVATURE OF MINIMUM SECTION 
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ONI LINE FOR CONSTANT VALUES OF RADIUS OF CURVATURE OF MINIMUM SECTION 


59 



1.950 


1.949 

1.948 
1.947 

1.946 
1.945 

1.944 
1.943 
1.942 
1.941 
1.940 
1.939 

1.938 
1.937 

.09 .10 .11 .12 .13 .14 .15 .16 .17 -18 

FIGURE A-9 VARIATION IN NON DIMENSION ALIZED MASS FLOW VERSUS 
DISTANCE BETWEEN MINIMUM SECTION AND MACH ONE LINE 
FOR VARIOUS CONSTANT VALUES OF THE RADIUS OF 
CURVATURE OF THE MINIMUM SECTION 



60 



BIBLIOGRAPHY 


1. McCarty, J. CCMD/HO, Letter of Transmittal, April 23, 1962, 
To: D. Thompson, M-P&VE-PA, Subject: Internal Charac- 
teristics. 


2. Oswatitsch, K. and Rothsteiii, W. , Flow Pattern in Converging- 
Diverging Nozzle, NACA-1215 

3. Sauer, R. , General Characteristics of the Flow Through Nozzles 
at Near Critical Speeds , NACA-TM-1147 

4. Shapiro, A. H. , The Dynamics and Thermodynamics of 
Compressible Fluid Flow, Volume II, The Ronald Press Company, 
1954 

5. Mendelson, R. S. , Method for Analysis of the Transonic Flow 
Field in the Throat of Axially Symmetric Rocket Nozzles, 
CCSD/HO TM-1069, June 1962 

6. Williams, J. R. , A Comprehensive Report on the Computer 
Program for the Design of a Converging-Diverging Rocket 
Nozzle, Allison Report RN 61-95 


61 


July 6, 1964 


APPROVAL 


TM X- 53084 


INTERIM REPORT ON METHODS OF DETERMINING 
THE TRANSONIC FLOW FIELD IN AN AXIALLY 
SYMMETRIC ROCKET NOZZLE 

By R. S. Mendelson 

The information in this report has been reviewed for security 
classification. Review of any information concerning Department 
of Defense or Atomic Energy Commission programs has been made 
by the MSFC Security Classification Officer. This report, in its 
entirety, has been determined to be unclassified. 



F. B. CLINE 


Acting Director, Propulsion and Vehicle Engineering Laboratory 


62 
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ABSTRACT 

Starting data are provided for internal characteristic programs 
pertaining to axially symmetric rocket nozzles. Several approximate 
methods are examined before a complete analysis is attempted. 

Using perfect frozen gas considerations, the potential function in 
the region of interest is approximated by a double power series in the 
space variables. The potential equation of motion, with the inviscid 
boundary condition, produces non-linear simultaneous equations. 

The non-linear equations are handled uniquely, and the results 
are utilized to describe the flow field. Different methods of checking 
the validity of the results are applied, and comparisons are made 
with other analyses. 

The remaining difficulties in the development of a production 
program that can handle arbitrary minimum section geometry and 
slightly varying mass flows are discussed. 
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